classdef SimonsQuasi  < Simons
    properties       
    end
    methods       
          
        function [Solution obj] = Run(obj)
            %% Culot Algorithm ;
            %% A quasi-Newton algorithm for first-orderSaddle Points
            %% PDF: -> 'D:\Saddle Points Literature\Transition States\Transition
            %% State\Deterministic Methods\Single ended\A quasi-Newton algorithm for first-orderSaddle Points.pdf'

            obj.Dims = 2;
            global X
%             X = [-1.0; 0.0]; %% zeros(1,Dims);
            %X = [0; 0];
            Solution = [];
            %X = [-1.5; -0.5];
            %% 1) Calculate the energy and the derivatives at the initial geometry
            [obj.E, obj.Grad, obj.H]= obj.CalculateRealEnergyDerivatives(X);
            isSolution = false;
            obj.Rk = 0.2;
            
            %% 2) Diagnolize the hessian matrix and choose an eigenvector
            %% vj
             if isnan(det(obj.H)) || isinf(det(obj.H)) || det(obj.H) ==0
                    fprintf('\n Singular Matrix encountered during the search \n');
                    return;
                end
            [obj.eigVectors eigValsM] = eig(obj.H);
            
            obj.eigVals = zeros(1,obj.Dims);
            for i = 1: obj.Dims
               obj.eigVals(i) =  eigValsM(i,i);
            end                
          
            obj.counter = 0 ;          
            for K = 1: obj.MaxIterations
                %% 3) Calculate Step Vector 
                [obj.n_pos, obj.n_zero, obj.n_neg] = obj.DetermineInertiaOfHessian(...
                    obj.eigVals, obj.Dims);
                obj.counter = obj.counter + 1;
                fprintf(' \n Current Step No : %d \t Solution: (%d,%d) \t ||g|| = %d #negEigenValues = %d \n', ...
                    obj.counter, X, norm(obj.Grad), obj.n_neg);
                Step = obj.MoveAStep();
                
                %% 5) Calculate the quadratic energy function approximation    
                %% 6) Calculate the energy at the point X + D upadte the trust radius
                E1 = obj.CalculateQuadraticEnergy(Step, obj.E, obj.Grad, obj.H);    
                
                Ce = obj.CalculatePredictedChangeInEnergy( obj.Grad, Step, obj.H);
                TrstRedius = TrustRegion;
                obj.Rk = TrstRedius.UpdateTrustRaduis(obj.Rk, E1, ...
                    obj.E, Ce, Step, obj.Grad, obj.H);                             
              
              
                %% is the current configruation is transition states or first order
                %% saddle points
                %% 7.2.1) Yes: report result
                
                GradRMS = norm(obj.Grad);
                MaxGrad = max(obj.Grad);
                MaxStep = max(Step);
                StepRMS = norm(Step);
                
                if StepRMS == Inf
                    continue;
                end

                isSolution = obj.isATransitionState(GradRMS, MaxGrad,...
                    MaxStep, StepRMS);
                if isSolution                     
                    break;
                end

                %% 7.2.2) No:
                %% Update the geometry X, Calculate the gradient vector at the new
                %% geometry, update the hessian matrix     
                X = X + Step;
                obj.E = E1;   
                Gk = obj.Grad;
                obj.Grad = Gradient(X);
                Gk1 = obj.Grad;
                obj.H = obj.UpdateHessian(Step, Gk1, Gk, obj.H);                                
                if isnan(det(obj.H)) || isinf(det(obj.H)) || det(obj.H) ==0
                    fprintf('\n Singular Matrix encountered during the search\n');
                    return;
                end
              
                [obj.eigVectors eigValsM] = eig(obj.H);         
                obj.eigVals = zeros(1,obj.Dims);
                for i = 1: obj.Dims
                   obj.eigVals(i) =  eigValsM(i,i);
                end                
                
            end
            if isSolution
                Solution = X;          
            end
            
        end        
        
        function Step = MoveAStep(obj)  
            Step = obj.CalculateQausiStep();
            if  obj.n_neg~= 1  || norm(Step) > obj.Rk 
                fprintf('\n Augemented QuasiNetwon Step ');               
                Step = obj.CalculateStep();
            end
            
        end
       
      
    end
end
    


